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An iterative algorithm for reconstructing the photon distribution from the random 
phase homodyne statistics is discussed. This method, derived from the maximum- 
likelihood approach, yields a positive definite estimate for the photon distribution 
with bounded statistical errors even if numerical compensation for the detector 
imperfection is applied. 

A fascinating topic studied extensively in quantum optics over past several years 
is the measurement of the quantum state of simple physical systems [1]. The central 
question is how to reconstruct a family of observables characterizing the quantum state 
from data that can be obtained using a feasible experimental scheme [2] . Most of the 
initial work on this topic was based on considerations of ideal, noise-free distributions 
of quantum observables that can be obtained only in the limit of an infinite number 
of measurements. However, any realistic experiment yields only a finite sample of 
data. Recognition of the full importance of this fact has led to the development of 
reconstruction techniques specifically designed to deal with finite ensembles of physical 
systems [3-6]. The main motivation for these developments is to optimize the amount 
of information on the quantum state that can be gained from a realistic measurement, 
and to distinguish clearly between the actual data obtained from an experiment and a 
priori assumptions used in the reconstruction scheme. 

One of ideas that have proved to be fruitful in the field of quantum state mea- 
surement is the principle of maximum-likelihood estimation. In particular, it has been 
recently applied to the reconstruction of the photon distribution from the random phase 
homodyne statistics [7] . The essential advantage of the maximum- likelihood technique 
is that physical constraints on the quantities to be determined can be consistently built 
into the reconstruction scheme. This reduces the overall statistical error, and auto- 
matically suppresses unphysical artifacts generated by standard linear reconstruction 
techniques [8], such as negative occupation probabilities of Fock states. However, this is 
achieved at a significantly higher numerical cost than that required by linear techniques. 
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Before we pass to the detailed discussion of the maximum-likelihood method, let 
us demonstrate its capability to improve the accuracy of the reconstruction. Fig. 1 
depicts the final result of processing Monte Carlo simulated homodyne data for a co- 
herent state and a squeezed vacuum state; both states have the average photon number 
equal to one. We have assumed imperfect photodetection characterized by a quantum 
efficiency rj — 85%, which is numerically compensated in the reconstruction scheme. 
It is seen that the maximum-likelihood estimate, marked with dots, is much closer to 
the exact distribution than probabilities obtained using the standard linear method of 
pattern functions. Moreover, purely artificial nonzero values for large photon numbers 
completely disappear when the maximum-likelihood method is applied. 

So, how does the maximum-likelihood method work, and what algorithm provides 
the bounded, positive definite estimate for the photon distribution? Let us start the 
discussion from considering the raw data recorded in an experiment. A single experi- 
mental run of a random phase homodyne setup yields a certain value of the quadrature 
observable q. This data, after an appropriate discretization, is recorded as an event in 
a vth bin. The probability of this event p v depends linearly on the photon distribution 
{£>«}: 

Pu{{Qn})=^A vn Q n . (1) 

n 

For a fixed n, the set of coefficients A vn describes the homodyne statistics for the nth 
Fock state. 

Repeating the measurement N times yields a frequency histogram {k n } specifying 
in how many runs the outcome was found in a ^th bin. These incomplete, finite data 
are the source of information on the photon distribution. The question is, how this 
information can be retrieved. The answer given by the maximum-likelihood estimation 
method is to pick up the photon distribution for which it was the most likely to obtain 
the actual result of the performed series of measurements. Mathematically, this is done 
by the maximization of the log- likelihood function [7] : 

C{{Q n }) = Y,ku^Pu{{Qn})-NY,Qn, (2) 

v n 

where N = ^„ k v is the total number of experimental runs. In the above formula, the 
method of Lagrange multipliers has been used to satisfy the condition that the sum of all 
probabilities is equal to one. As the estimate for {g n } is supposed to describe a possible 
physical situation, the search for the maximum likelihood is o priori restricted to the 
manifold of real probability distributions. Geometrically, this manifold is a simplex 
defined by the set of inequalities: 

g n >0, n = 0,1,2,... 

8n = 1- 

n 

Thus, the physical constraints on the reconstructed quantities are naturally incorpo- 
rated in the maximum-likelihood reconstruction scheme. 
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Fig. 1. Reconstruction of the photon distribution from Monte Carlo homodyne data. The bars 
depict exact photon distributions for a coherent state (top) and a squeezed vacuum state (bot- 
tom), both the states with the average photon number equal to one. For each of these states, 
N = 10 5 homodyne events were simulated. The events were divided into 100 bins covering the 
interval —5 < q < 5. The simulated data were used to reconstruct the photon distributions 
using the maximum-likelihood method (•) and the linear pattern function technique ({>). In 
the latter case, some of reconstructed g n 's are beyond the scale of the graphs. Note that for 
small n, Fock state occupation probabilities reconstructed using both the methods are very 
close, and the symbols • and <0 overlap. The maximum-likelihood estimates were obtained 
from 8000 iterations of the EM algorithm, starting from a uniform distribution for < n < 20. 
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The reconstruction of the photon distribution formulated in the maximum-likelihood 
approach belongs a very wide class of linear inverse problems with positivity con- 
straints [9] . This class encompasses a variety of problems appearing in as diverse fields as 
medical imaging and financial markets. As discussed by Vardi and Lee [9], a straight- 
forward method for solving these problems is provided by the so-called expectation- 
maximization (EM) algorithm. In the following, we will present a heuristic derivation 
of this algorithm applied the reconstruction of the photon distribution. 

Let us consider the necessary condition for the maximum of the function C({g n }). 
The partial derivatives of the log-likelihood function are given by: 

dQm ^ Pv{{Qn\) 

For each m, the partial derivative dC/dg m must vanish unless the maximum is located 
on a face of the simplex for which g m = 0. Thus we have that either dC/dg m = or 
g m = 0. These two possibilities can be written jointly as 



{{Qn}) 



— AM = 0, m = 0,1,2,.... (4) 



It is convenient to rearrange this condition to the form which defines the maximum- 
likelihood estimate as a fixed point of a certain nonlinear transformation. Such a form 
immediately suggests an iterative procedure for finding the fixed point by a multiple 
application of the derived transformation: 

g& +1) = gffE77 ^TS ■ m = 0,1,2,... (5) 

where upper indices in parentheses denote consecutive approximations for the photon 
distribution. 

In fact, Eq. (5) provides a ready-to-use iterative method for reconstructing the 
photon distribution, which is a special case of the EM algorithm [9] . Given an approxi- 
mation of the photon distribution {g$}, the next one is simply evaluated according to 
the right hand side of Eq. (5). When sufficient mathematical conditions arc fulfilled, 
this procedure converges to the maximum-likelihood solution. Thus, the complex prob- 
lem of constrained multidimensional optimization is effectively solved with the help of a 
simple iterative algorithm. The maximum-likelihood estimates depicted in Fig. 1 were 
obtained from 8000 iterations of this algorithm. The initial distributions were assumed 
to be uniform for photon numbers n < 20 and equal to zero above this cut-off value. 

An essential yet delicate matter in quantum state measurements is the role played by 
the detection efficiency. The impact of imperfect detection on the maximum-likelihood 
reconstruction scheme can be understood using a simple intuitive argument. Accord- 
ing to Eq. (1), the homodyne statistics is a sum of components generated by the Fock 
states |n), with the weights given by the appropriate occupation probabilities g n . The 
shape of each component is described by the set of coefficients A un , with n treated as 
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Fig. 2. Components of the random phase homodyne statistics generated by different Fock 
states, for decreasing efficiency r\ of the homodyne detector. The homodyne statistics of an 
arbitrary quantum state is a sum of these components with weights defined by the photon 
distribution {g n }- 

a fixed parameter. Fig. 2 depicts several of these components for different values of the 
detector efficiency rj. In the case of unit efficiency, the contribution originating from 
the Fock state \n) is given by the squared modulus of the nth energy eigenfunction 
in the position representation, and exhibits characteristic oscillations. The important 
point is that each of these contributions has a unique location of maxima and minima. 
Thus, each number state leaves a specific, easily distinguishable trace in the homodyne 
statistics. Roughly speaking, this is what allows the maximum-likelihood method to 
resolve clearly the contributions from different number states. When the detection is 
imperfect, homodyne statistics generated by Fock states become blurred, and oscilla- 
tions disappear. This makes the shape of the contributions from higher Fock states quite 
similar. Consequently, it becomes more difficult for the maximum- likelihood method 
to resolve components generated by different number states. Of course, this intuitive 
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argument should be supported by quantitative considerations. Mathematically, the ef- 
fect of imperfect detection is reflected by the shape of the log-likelihood function C. 
For low detection efficiency, it becomes flatter, and its maximum is not sharply peaked. 
This increases the statistical uncertainty of the reconstructed photon distribution, and 
results in a slower convergence of the EM algorithm. 

Finally, let us note that the maximum-likelihood algorithm docs not make any use 
of the specific form of the coefficients A vn linking the photon distribution with the 
homodyne statistics. In fact, it can be applied to any phase-insensitive measurement 
whose statistical outcome depends on the photon distribution via a linear relation of 
the form assumed in Eq. (1). 
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